########################
# Robustness Check 
########################
library(tidyverse)
library(stargazer)
library(ggeffects)
library(emmeans)

# Set your working directory
setwd("~/Dropbox/2001_replication/draft/data")

# Load data
load("output/final_results.RData")
dat <- read.csv("output/rueda_rep.csv")

# Make sure that year and occupation are factor variables
dat$year <- as.factor(dat$year)
dat$euroesec <- as.factor(dat$euroesec)

#  1. Logistic regression with country fixed effects

# Southern European Countries
fit_southern_fe <- glm(redpref1 ~ incomedif*forpop + age + gender +
                        education + religion + union + euroesec +
                        socgdp + year + cntry,
                      data = dat %>% filter(brncntr == 1, welfare == "southern"),
                      family = binomial(link = logit))
summary(fit_southern_fe)

# Non-Southern European Countries
fit_nonsouthern_fe <- glm(redpref1 ~ incomedif*forpop + age + gender +
                            education + religion + union + euroesec +
                            socgdp + year + cntry,
                          data = dat %>% filter(brncntr == 1, welfare != "southern"),
                          family = binomial(link = logit))
summary(fit_nonsouthern_fe)

# Nordic countries (social democratic)
fit_socdem_fe <- glm(redpref1 ~ incomedif*forpop + age + gender +
                       education + religion + union + euroesec +
                       socgdp + year + cntry,
                     data = dat %>% filter(brncntr == 1, welfare == "social_dem"),
                     family = binomial(link = logit))
summary(fit_socdem_fe)

# Continental European countries (corporatist)
fit_corp_fe <- glm(redpref1 ~ incomedif*forpop + age + gender +
                  education + religion + union + euroesec +
                  socgdp + year + cntry,
                data = dat %>% filter(brncntr == 1, welfare == "corporatist"),
                family = binomial(link = logit))
summary(fit_corp_fe)

# Liberal countries
fit_lib_fe <- glm(redpref1 ~ incomedif*forpop + age + gender +
                 education + religion + union + euroesec +
                 socgdp + year + cntry,
               data = dat %>% filter(brncntr == 1, welfare == "liberal"),
               family = binomial(link = logit))
summary(fit_lib_fe)

# Table 2 in the Appendix
stargazer(fit_southern_fe, fit_nonsouthern_fe,
          fit_socdem_fe, fit_corp_fe, fit_lib_fe, type = "latex")

## 2. Predicted probabilities averaged over different levels of categorical variables
# See: https://strengejacke.github.io/ggeffects/articles/technical_differencepredictemmeans.html 

pred_southern <- ggemmeans(fit_southern, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_southern)

pred_socdem <- ggemmeans(fit_socdem, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_socdem)

pred_corp <- ggemmeans(fit_corp, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_corp)

pred_lib <- ggemmeans(fit_lib, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_lib)

pred_southern$facet <- "Southern"
pred_socdem$facet <- "Nordic"
pred_corp$facet <- "Continental"
pred_lib$facet <- "Liberal"
pred_regime <- rbind(pred_southern, pred_socdem, pred_corp, pred_lib)

# Figure 3
pred_plot(pred_regime)

pred_mcwr <- ggemmeans(fit_d_mcwr, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr"))
pred_mcwr$facet <- ifelse(pred_mcwr$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")

# Figure 4
pred_plot(pred_mcwr)

# Figure 5
pred_childcare <- ggemmeans(fit_childcare, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr")) 
pred_childcare$facet <- ifelse(pred_childcare$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_childcare)

# Figure 6
pred_unemp <- ggemmeans(fit_unemp, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr"))
pred_unemp$facet <- ifelse(pred_unemp$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_unemp)

# Figure 7
pred_imm1 <- ggemmeans(fit_imm1, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr")) 
pred_imm1$facet <- ifelse(pred_imm1$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_imm1, ylab = "Support for Immigration")

# Figure 8
pred_imm2 <- ggemmeans(fit_imm2, c("incomedif [-2.5:9.2]", "forpop[5.3, 14.7]", "d_mcwr")) 
pred_imm2$facet <- ifelse(pred_imm2$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_imm2, ylab = "Support for Immigration")
